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ABSTRACT 

We investigate the effect of including diffuse field radiation when modelling the radia- 
tively driven implosion of a Bonnor-Ebert sphere (BES). Radiation-hydrodynamical 
calculations are performed by using operator splitting to combine Monte Carlo pho- 
toionization with grid-based Eulerian hydrodynamics that includes self-gravity. It is 
found that the diffuse field has a significant effect on the nature of radiatively driven 
collapse which is strongly coupled to the strength of the driving shock that is estab- 
lished before impacting the BES. This can result in either slower or more rapid star 
formation than expected using the on-the-spot approximation depending on the dis- 
tance of the BES from the source object. As well as directly compressing the BES, 
stronger shocks increase the thickness and density in the shell of accumulated material, 
which leads to short, strong, photo-evaporative ejections that reinforce the compres- 
sion whenever it slows. This happens particularly effectively when the diffuse field is 
included as rocket motion is induced over a larger area of the shell surface. The forma- 
tion and evolution of 'elephant trunks' via instability is also found to vary significantly 
when the diffuse field is included. Since the perturbations that seed instabilities are 
smeared out elephant trunks form less readily and, once formed, are exposed to en- 
hanced thermal compression. 

Key words: stars: formation - ISM: HII regions - ISM: kinematics and dynamics - 
hydrodynamics - radiative transfer - methods: numerical 



1 INTRODUCTION 

The majority of stars form in clusters, situated in molecular 
clouds that range in size from less than a single parsec to 
several hundred parsecs (Lada Lada 2003 ). In order for 
star formation to occur, gravitational collapse of material 
has to overcome internal thermal pressure, supersonic 
material motions (turbulence) and magnetic fields (see e.g. 
[Preibisch Zinnecke r|2007||Hartmannf 2QQ9). The presence 
of OB stars in these systems has a dramatic impact on 
the surrounding material (and therefore star formation), as 
they emit large amounts of high energy radiation that pho- 
toionizes gas and gives rise to propagating ionization and 
shock fronts ( Elmegr een|20 11 ). They also inject mechanical 
energy into the surroundings in the form of high-speed 
stellar winds and, eventually, supernova explosions. The 
net impact of radiative feedback from massive stars on 
star formation efficiency in a molecular cloud is currently 
unclear, though a number of individual processes that either 
inhibit or induce further star formation has been identified. 
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The two main radiative feedback mechanisms that 
inhibit star formation are the dispersal of material that 
might otherwise move towards the centre of the molecular 



clouds' gravitational potential (e.g. Herbig 1962) and the 
possibility of driving and maintaining turbulence that sup- 
ports against collapse (e.g. Peters et al.||2008" Gritschneder 
|etral.||2Q09t . 

The two primary established mechanisms for the in- 
duction of star formation are consequences of the expanding 
ionization and shock fronts about a massive star. The first 
is 'collect and collapse' (Elmegreen k, Lada|1977 Elmegreen 



et al.|1995 Dale et al.|2007 ) in which material that is accu 
mulated by the expanding ionization front of an Hll-region 
becomes locally gravitationally unstable, fragmenting and 
collapsing to form stars. This mechanism is supported ob- 
servationally by, for example, the identification of massive 
fragments located in a dust ring surrounding the HII region 
ROW 79 bylZavagno et al.l (12006). 



The second mechanism is radiatively driven implosion 
(RDI), ( Bertoldi]|1989 ) in which radiatively induced shocks 
drive into otherwise stable pre-existing density structures 
and cause them to collapse and form stars. RDI has been 
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modelled by various groups, for example |Kessel-Deynet 


& Burkert ( 


2003), Arthur k Hoare (2006), |Gritsclineder 


et al. (2009 


), Henney et al. (2009), Bisbas et al. (2011) 



and produces objects similar to observed bright rimmed 
clouds (BRCs) in HII regions ( |Ogura et aL]|2002t . The rate 
at which collapse occurs and the associated star forming 
efficiency of these RDI models has unsurprisingly been 
found to be very sensitive to the incoming flux. 

A variant of collect and collapse, in which the trigger 
for fragmentation of an ionization front is one of the 
many possible hydrodynamic or radiation hydrodynamic 
instabilities has also been explored (e.g. 
uarcia- Seg ura Franco] [19961 IWilhamsl 
et al. 2006). When sufficiently perturbed, faster moving 
components of the ionization front will funnel material 
transversely to the direction of I-front propagation, deposit- 



Vishniacl 



|2002| [Mizutal 



1983 



ing it in the path of slower moving components (Vishniac 



1983 Garcia-Segura & Franco 1996). This results in a 



collection of pillar-like objects with dense tips, much as 
is observed in e.g. the pillars of creation in the Eagle 
Nebula or the Dancing Queen's Trunk in NGC7822 (see e.g. 
Schneps et al.|1980| [Reach et al.|20"04||Chauhan et al.|2011 ). 



Recently, [Gritschneder et al. ( 2009 ) and Gritschneder 
et al. ( 2010| used the ray-tracing, smooth particle hydro- 



dynamics (SPH) code, iVINE ( [Gritschneder et al.|[2009t to 
consider the effects of a radiation front impinging upon a 
turbulent neutral medium. They found that radiation could 
support turbulence, which would prohibit large-scale star 
formation by supporting against collapse. They also found 
that ionizing radiation rapidly penetrated lower density 
regions, heating them up and compressing the remaining 
higher density structures. This again resulted in pillars 
with dense cores at the tips where stars might form. They 
termed this process 'radiative round-up'. 

Each of the aforementioned processes has been the sub- 
ject of numerical modelling. Due to the intensive computa- 
tional demand of these models, at least in three dimensions, 
a number of approximations have necessarily been devel- 
oped. [ErcoIano^GTitsch^ (2011a) provides an overview 
and evaluation of some of the main approximations, sum- 
marised as follows. 

a) Considering a monochromatic radiation field, usu- 
ally Lyman 13.6eV photons, reduces calculation timescales 
(as with all of these approximations) compared with poly- 
chromatic models. This is because the entire source spec- 
trum does not need to be resolved and a single value can 
be used for the gas opacity. However, this speed up is at 
the expense of being able to reliably calculate the resulting 
ionization and temperature structure of the system and ne- 
glects effects due to radiation hardening (though radiation 
hardening effects can be estimated in monochromatic codes, 
e.g. |Mellema et al.|[2006| . 

b) Use of simplified thermal balance calculations, for 
example calculating the temperature as a simple function 
of the ionization fraction. This is more straightforward to 
implement and results in faster calculations than solving the 
thermal balance by comparing heating and cooling rates. 

c) Assuming that the system is in photoionization equi- 



librium. This is valid where recombination timescales are 
shorter than the dynamical timescales of the gas. 

d) The 'on the spot' (OTS) approximation. Under this 
scheme diffuse field photons, those generated in recombi- 
nation events, are not treated. This is justified in regimes 
where diffuse field photons will not propagate very far 
and therefore not modify the global ionization structure 
significantly. It is however, questionable in regions of low 
or rapidly varying density. Because diffuse field photons 
are emitted isotropically, it is possible that shadowed 
regions will not be exposed to a realistic amount of ioniz- 
ing radiation when modelled under the OTS approximation. 



It is not yet clear what impact these approximations 
have on both the inhibiting and inducing mechanisms 
described above. In an effort to understand the effect of the 
diffuse field, Ercolano Gritschneder ( [2011a ) compared 
snapshots from the radiative round-up models run by 
iVINE, with full radiative transfer calculations using the 
Monte-Carlo radiative transfer code MOCASSIN ( [Ercolano" 
et al. 2003^ and noted significant differences between the 



ionization and temperature structures calculated by the 



two codes. Ercolano & Gritschneder (2011b) subsequently 



attempted to account for the thermal effects of the diffuse 
field in iVINE by identifying shadowed regions and assign- 
ing them a parameterised temperature as a function of 
density based on comparisons of iVINE and MOCASSIN 
snapshots. Using this shadowing scheme a significant effect 
on the resulting pillar structures was observed. Far fewer 
pillars were formed and those that remained were narrower, 
denser and often cut off from the parent molecular cloud 
due to their paramaterised increased exposure to ionizing 
radiation from the diffuse field. This leads to earlier 
triggering of star formation and reduces the efficiency of 
radiation driven turbulence. This shadowing scheme is not 
without drawbacks, being susceptible to erroneous heating 
of true shadowed regions. The spread in temperatures 
at a given density calculated by MOCASSIN also results 
in quoted typical errors in parameterized temperature of 
approximately 50%. It is certainly clear from this work that 
a more comprehensive knowledge of the effects of the diffuse 
field would be valuable. It will be necessary to establish 
just how different the result of a radiation hydrodynamics 
calculation can be when incorporating the diffuse field di- 
rectly, to validate or reassess the use of simplified radiation 
handling in radiative feedback simulations. 



In this work we use the radiation hydrodynamics code 
TORUS ([Harries[[2000l [Kurosawa et al.[[2004| [Harries et al 



2004 Acreman et al. 2010 Harries 2011) to investigate 
the effects on radiatively driven implosion of a more so- 
phisticated treatment of the diffuse field than previously 
applied in a radiation hydrodynamics calculation. Specif- 
ically, we will systematically deduce the relative effects 
of using a monochromatic OTS, polychromatic OTS and 
polychromatic-diffuse radiation field on the overall nature 
of collapse. 
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2 NUMERICAL METHOD 

2.1 Hydrodynamics 

We use a flux conserving, finite volume hydrodynamics al- 
gorithm. It is total variation diminishing (TVD) and makes 
use of the Superbee flux limiter ( Roe|1985 ). We also incorpo- 
rate a Rhie-Chow interpolation scheme (|R hie ChowjlOSS ) 
to avoid 'checkerboard' effects that may otherwise appear in 
cell-centred grid-based hydrodynamics codes. 

2.2 Self-Gravity 

In order to include self- gravity we solve Poisson's equation 



(1) 



where (f) and p are the gravitational potential and matter 
density respectively. 

We use a multigrid method employing Gauss-Seidel 
sweeps with successive over-relaxation to solve a linearized 
version of equation ^ over the grid. We impose Dirichlet 
boundary conditions with the potential at the boundary set 
to a multipole expansion of the matter interior to the bound- 
ary (including terms up to the quadrupole). The gravita- 
tional potential is subsequently included as a source term in 
the hydro dynamical equations. 



2.3 Photoionization 

Photoionization calculations are performed here using an 
iterative Monte Carlo photon packet propagating routine, 



similar to that of Ercolano et al. (2003) and Wood et al 



2004| which in turn are based on the methods presented by 



Lucy (1999 



Photon packets are collections of photons for which the 
total energy e remains constant, but the number of pho- 
tons contained varies for different frequencies These are 
initiated at stars in the model, with frequencies selected ran- 
domly based on the emission spectrum of the star. The con- 
stant energy value e for each photon packet is simply the 
total energy emitted by stars (luminosity L) during the du- 
ration At of the iteration divided by the total number of 
photon packets N: 

LAt 



(2) 



Once initiated, a photon packet will propagate for a 
path length p determined by a randomly selected optical 
depth as detailed in Harries & Howarth| (|1997|) . If the pho- 
ton packet fails to escape a cell after travelling p then its 
propagation ceases and an absorption event occurs. Under 
the OTS approximation, once a photon packet is absorbed 
it is ignored, being assumed to either have been re-emitted 
with a frequency lower than that required for photoioniza- 
tion, or provide negligible further contribution to the ion- 
ization structure by causing further photoionization on only 
small scales. Using the principle of detailed balance, after an 
absorption event a new photon packet is immediately emit- 
ted from the same location with a new isotropically random 
direction and a new random frequency based on the tem- 
perature dependent emission spectrum. This process repeats 
until the photon packet escapes the grid, or its propagation 
time matches the current simulation time. For each cell in 



the grid the sum of the paths / that photon packets traverse 
is recorded. 

Note that the energy density dU of a radiation field is 
given by 



(3) 



where c is the speed of light. A photon packet having a path 
/ in a particular cell contributes an energy e(//c)/At to the 
time-averaged energy density of that cell. Thus by summing 
over all paths / the energy density of a given cell (volume 
V) can be determined, and by equation [s] the mean intensity 
may be estimated: 



47rJ^ 



cAtV 



(4) 



This is then used to obtain ionization fractions by solving 
the ionization balance equation ( Osterbrock||1989 ) 

n{X'+^) _ 1 



4:7TJuau{X'^)diy 



(5) 



where n{X^), a{X^), au{X^), rie and are the number den- 
sity of the i^^ ionization state of species X, recombination 
coefficient, absorption cross section, electron number den- 
sity and the threshold frequency for ionization of species X* 
respectively. In terms of Monte Carlo estimators (equation 
|4|, equation [5] is given by 

n{X'+^) _ e 



n{X^) AtVa{X^)n, 



E 



la^{X' 



(6) 



This approach has the advantage that photon packets con- 
tribute to the estimate of the radiation field without having 
to undergo absorption events, thus even very optically thin 
regions are properly sampled. 

TORUS performs photoionization calculations that in- 
corporate a range of atomic species and in which thermal 
balance in each cell is calculated by iterating on the tem- 
perature until the heating and cooling rates match. 

For the radiation hydrodynamics models in this paper a 
simplified thermal balance calculation is used and the only 
species considered are atomic and ionized hydrogen. This is 
to allow for comparison with previous works that use these 



schemes such as Gritschneder et al. ( 2009 ) and Bisbas et al. 
(2009). The temperature is calculated by interpolating be- 
tween pre-determined temperatures, Tn and Tio, ascribed to 
the state of fully neutral and fully ionized gas respectively as 
a function of the newly calculated fraction of ionized atomic 
hydrogen in the i^^ cell 77^: 



T, =Tn + ry^(Tio - Tn) 



(7) 



In this paper we use Tn = 10 K and Tio = 10000 K for all 
models that use this simplified thermal balance calculation. 



2.4 Radiation Hydrodynamics 

The hydrodynamics and photoionization schemes outlined in 
sections [2rT] and (23] are combined using operator splitting to 
perform radiation hydrodynamics calculations. A photoion- 
ization calculation is initially run to convergence, this gener- 
ally allows subsequent calculations to run relatively quickly 
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given that they are usually minor perturbations of the pre- 
vious state. We then perform photoionization and hydrody- 
namics steps sequentially. In this work, the photoionization 
calculation for a time step is always calculated prior to the 
hydrodynamics calculation. 

This operator splitting technique is flexible and rela- 
tively conceptually straightforward. It is however very com- 
putationally expensive, requiring a large number of Monte- 
Carlo photoionization calculations that render the gravita- 
tional and hydrodynamic components of the calculation neg- 
ligible in comparative computational cost. Fortunately, this 
approach can be efficiently parallelized. 



2.5 Implementation 

Despite the high computational cost of this radiation hy- 
drodynamics scheme, it is extremely scalable. The compu- 
tational domain is decomposed into subdomains over which 
the components of the radiation hydrodynamics calculation 
are computed by an individual processor (thread). In three 
dimensions subdomains take the form of cubes of equal vol- 
ume, which at present can be either 1/8, 1/64 or 1/512 of 
the total domain volume. An additional master thread per- 
forms governing and collating operations, giving a total of 65 
threads for each of the three-dimensional models performed 
here. At lower dimensionality the grid can be decomposed 
in a similar manner into equally sized squares (2D) or lines 
(ID). 

In the photoionization component of a calculation, pho- 
ton packets are communicated between threads in stacks 
rather than individually to reduce the communication la- 
tency overhead. TORUS stores quantities using an octree 
AMR grid however adaptive refinement and coarsening of 
the grid is not yet fully tested in the hydrodynamics rou- 
tine, a fixed grid is therefore used for models in this work. 

The RDI calculations presented in this paper were run 
on an SGI Altix ICE system using 65 2.83GHz Intel Xeon 
cores across 9 dual quad-core compute nodes. These typ- 
ically completed within 600-1000 hours of wall time. The 
models that include the diffuse field do not necessarily take 
the longest time to complete, as the additional hydrody- 
namic and photoionization calculations required for models 
that develop the highest velocity material motions outweigh 
the additional time taken for each photoionization calcula- 
tion when the diffuse field is included. 



3 NUMERICAL TESTS 

A number of tests have been conducted to confirm that 
TORUS is in agreement with accepted benchmarks. We 
treat the photoionization, hydrodynamics and self-gravity 
in isolation as well as the radiation hydrodynamics. Com- 
prehensive tables of test parameters are provided to enable 
replication by other codes. 



3.1 Photoionization: The HII40 Lexington 
Benchmark 

The HII40 Lexington benchmark is a one dimensional test in 
which the equilibrium temperature and ionization structure 
of an HII region heated by a star at 40000 K is calculated 



Table 1. Lexington benchmark parameters. 



Variable (Unit) 


Value 


Description 


Teff(K) 


40000 


Source effective temperature 




18.67 


Source radius 


nil (cm-3) 


100 


Hydrogen number density 


logio(He/H) 


-1 


Helium abundance 


logio(C/H) 


-3.66 


Carbon abundance 


logio(N/H) 


-4.40 


Nitrogen abundance 


logio(0/H) 


-3.48 


Oxygen abundance 


logio(Ne/H) 


-4.30 


Neon abundance 


logio(S/H) 


-5.05 


Sulphur abundance 


L (cm) 


4.4x1019 


Computational domain size 


'^cells 


1024 


Number of grid cells 



and compared with the output of one of the many codes 
that reproduce the accepted result (see Ferland|1995 ). Here 
we calculate a comparison set of results using the one di- 
mensional semi-analytic code Cloudy ( Ferland et al.|p^98 ), 
one of the original contributors to the benchmark. The sys- 
tem modelled comprises a star at 40000 K at the left hand 
edge of the grid, size 4.4xl0^^cm comprising 1024 cells. This 
test incorporates more species than only hydrogen and does 
not rely on the simplified thermal balance calculation of 
equation [7| rather the temperature of the cell is determined 
through comparison of the heating and cooling rates. It also 
includes treatment of the diffuse field. 

A full list of parameters used for this benchmark is given 
in Table [l] and the resulting temperature and ionization 
fractions as calculated using both TORUS and Cloudy are 
shown in Figure [l] TORUS is consistent with the Cloudy 
temperature distribution to within 10% and is generally 
much better than this. The higher temperature calculated 
by TORUS in the inner regions is in agreement with the 
result obtained in Wood et al. (2004). The hydrogen and 
helium ionization fractions agree extremely well, this is of 
particular importance with regard to the hydrogen-only ra- 
diation hydrodynamics models in this paper. The other ions 
match to within similar levels of agreement as |Ercolano et aL\ 
( |2003| and [Wood et al. ( |2004| . Discrepancies in the result of 
this benchmark are usually attributed to differences in the 
atomic data used by the codes that are being compared. 



3.2 Hydrodynamics Tests 

3.2.1 Sod Shock Tube 

The Sod shock tube test is a simple 1-dimensional model, 
initially comprising two equal volumes separated by a par- 
tition. Both partitions, the left hand partition (LHP) and 
right hand partition (RHP), contain ideal gases at zero ve- 
locity with different densities and hence different initial pres- 
sure and energy. At time t = Os, the partition is removed 
and the system allowed to evolve. The state after a time t 



is detailed by |Sod| ( [1978 ). 

We use an adiabatic equation of state, with adiabatic 
index 7 = 7/5. Unless otherwise stated, a value of 0.3 is used 
for the Courant-Friedrichs-Lewy (CFL) parameter in all hy- 
drodynamics models. In this model boundary conditions are 
reflective and 1024 cells are used for the grid. 

A summary of the model parameters is given in Ta- 
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Figure 1. Top Left: The temperature distribution for the Lexington benchmark. Top Right: Hydrogen and hehum ionization fractions. 
Middle Left: Oxygen ionization fractions. Middle Right: Carbon ionization fractions. Bottom Left: Nitrogen ionization fractions. Bottom 
Right: Neon ionization fractions. 



ble |2] and the result as computed by TORUS at t = 0.2 s 
is shown in Figure |2] The features in this result are, from 
right to left; the initial density in the RHP, a shock wave 
which forms as a result of the low density material recoiling 
away from the high density material, a contact discontinu- 
ity between the high density material of the LHP and the 
low density material of the RHP, a rarefaction wave formed 



because the contact discontinuity acts as a piston drawing 
left hand material to the right and the initial density of the 
LHP. 

TORUS shows excellent agreement with the analyti- 
cal solution. The slight density dip at the rarefaction wave- 
contact discontinuity interface arises because TVD flux lim- 
ited schemes only smooth out oscillations near sharp shocks. 
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Table 2. Sod shock tube parameters. 



Variable (Unit) 


Value 


Description 






0.8 


pi (gcm-3) 


1 


Initial density in LHP 








P2 (gcm-3) 


0.125 


Initial density in RHP 








7 


7/5 


Adiabatic index 




m 

'E 


0.6 


Pi (dyncm~^) 


1 


Initial pressure in LHP 




u 




P2 (dyncm~^) 


0.1 


Initial pressure in RHP 








El (ergcm"'^) 


2.5 


Initial energy density in 


LHP 


C 
(U 

Q 


0.4 


El (ergcm""^) 


2.0 


Initial energy density in 


RHP 






E.O.S. 


Adiabatic 


Equation of state 








L (cm) 


1 


Computational domain 


size 




0.2 


'^cells 


1024 


Number of grid cells 
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Table 3. Sedov- Taylor parameters. 



Variable (Unit) 


Value 


Description 


7 


7/5 


Adiabatic index 


V (cms"-*-) 





Initial velocity 


p (gcm-2) 


1 


Initial surface density 


n (cm) 


0.01 


Energy dump zone radius 


EBlast (ergcm-2) 


3183.1 


Dump zone surface energy density 


Eo (ergcm~2) 


1 X 10-^ 


Ambient surface energy density 


E.O.S. 


Adiabatic 


Equation of state 


L(cm) 


1 


Computational domain size 


''^cells 


5122 


Number of grid cells 



0.0 0.2 0.4 0.6 0.8 l.( 

Distance, cm 

Figure 2. Density distribution of the Sod shock tube test at 
t = 0.2s, showing both the result given by TORUS (blue crosses) 
and the analytical solution (red line). 



Taylor] |l950bj and .Sedo^ p946 ). 

The initial ratio of thermal energy in the circular region 
to the rest of the grid is 3 x 10^ : 1. Given the extreme nature 
of this model, a CFL parameter value of 0.08 was required 
in order to capture the early stages of evolution without 
numerical instability arising. The boundary conditions used 
in this test are all reflective and the grid comprises 512^ 
cells. A full table of parameters is given in Table [3] and a 
comparison of the density distribution at time t = 0.03 s, as 
calculated both by TORUS and analytically, is given in Fig- 
ure |3] TORUS demonstrates a good level of agreement with 
the analytical solution but, as with all numerical schemes, 
suffers from numerical diffusion. This is responsible for the 
slight broadening of the shock and reduction in the peak 
amplitude compared to the analytical solution. 



3.2.3 Kelvin- Helmholtz Instability 

Kelvin-Helmholtz (KH) instabilities are vortices that form 
at an interface between two materials due to shear forces 
(Von Helmholtz||1868| |Kelvin||1871|). Adaptations to SPH 



This is a necessary compromise, so that existing physical os- 
cillations are not unphysically damped. 



3.2.2 Sedov- Taylor Blast Wave 

The Sedov- Taylor blast wave is an extreme model, which 
tests the advection scheme beyond the demands that will 
be made of it in star formation applications. In this 2- 
dimensional test a large amount of energy is injected into 
a circular region, radius 0.01 cm, of a constant density ideal 
gas, causing a blast wave. Self-similar analytical solutions 
for the time evolution of such a blast wave were found by 




TORUS 
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Figure 3. Density distribution of the Sedov- Taylor blast wave at 
t = 0.03 s showing both the result given by TORUS (blue crosses) 
and the analytical solution (red line). 



codes have recently found to be required in order to form 
KH instabilities by [Agertz et ah] ( |2007| ) and |Price| ( |2008| ), 
thus reproducing these has proved to be an important test 
of hydrodynamical algorithms. 

The 2-dimensonal system modelled here follows |Price| 
(2008|, comprising two fluids in contact at different density 



and velocity, such that the ratio of their densities is 2:1 and 
their velocities are equal in magnitude but in opposite di- 
rections. Periodic boundary conditions are used at the ±x 
boundaries and reflective conditions at the ±z boundaries, 
corresponding to housing the system in a pipe that extends 
indefinitely in dbx. At time t = Os the interface between 
these fluids is subject to a perturbation of the form 
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Table 4. Kelvin-Helmholtz parameters. 



Variable (Unit) 


Value 


Description 


7 


5/3 


Adiabatic Index 


Pi (gcm-2) 


1 


Ambient fluid initial surface density 


P2 (gcm-2) 


2 


Central fluid initial surface density 


ui (cms"-*-) 


-0.5 


Ambient fluid initial velocity 


U2 (cms"-*-) 


+0.5 


Central fluid initial velocity 


E.O.S 


Adiabatic 


Equation of State 


A 


0.025 


Constant in perturbation equation 


A (cm) 


1/6 


Wavelength of perturbation 


L (cm) 


1 


Computational domain size 


'^cells 


5122 


Number of grid cells 




Figure 4. Surface density distribution in gem 2 at tkh across 
the lcm2 domain of the Kelvin-Helmholtz instability model. 



r Asin(-27r(x+ l/2)(l/6)), |z - 0.25| < 0.025 
^ - <^ (8) 
I Asin(27r(x + l/2)(l/6)), \z + 0.25| < 0.025. 



vortices should then 
timescale given by 

27r 

Tkh = — . 



form within a characteristic KH 



(9) 



Where, for materials in contact with density pi and p2 and 
velocities vi and V2 subject to a periodic perturbation of 
wavelength A 



27r {pip2 



A 12 



V\-V2\ 



A 



(10) 



A full table of parameters used for this test is given in Ta- 
ble ^ Using these parameters and equations [9] and ^] we 
obtain a KH timescale of approximately 0.35 seconds, the 
time within which vortices should form. A plot of the den- 
sity distribution as calculated by TORUS at tkh is given 
in Figure [4] and clearly demonstrates that primary and sec- 
ondary vortices have formed within the KH timescale. This 
test has also been successfully performed using density ratios 
of 5:1 and 10:1. 




Figure 5. Surface density distribution in gcm~2 showing the 
'mushroom' formed in our 1 cm2 domain size, Ray leigh- Taylor 
instability model. It comprises a penetrating column of material 
with Kelvin-Helmholtz instabilities forming where the materials 
flow past one another. 



Rayleigh- Taylor Instability 

Ray leigh- Taylor instabilities can arise where a material is 
'on top' of a second lower density material in a gravitational 
potential and the interface between them is subject to a 
perturbation fRay leigh 1900: Tay lor] 1950a '). They are man- 
ifested as 'Ray leigh- Taylor fingers' that propagate into the 
low density material along which Kelvin-Helmholtz instabil- 
ities may form. This gives rise to a characteristic mushroom 
shape. In the following test we have selected system param- 
eters that should give rise to this characteristic structure. 

At time t = s the interface between two different den- 
sity materials in a 1 cm^ box in the presence of a gravita- 
tional field is subject to a small disturbance of magnitude 
—0.055 cm s~^ across a finite range 0.45 < x < 0.55 of the 
interface, the system is then left to evolve. We use periodic 
boundary conditions at the =bx-direction bounds and reflec- 
tive at the iz-direction bounds. The gravitational potential 
(j) at height z is given by 



(z) gz 



(11) 



where g is the acceleration due to gravity, here equal to 
0.1 cm s~^. 

A summary of parameters used in this test is given in 
Table [5] Figure |5] shows the distinctive mushroom-shape 
formed via this method at time t = 5s. The main body 
of the mushroom is the Rayleigh- Taylor finger, at the tip of 
which Kelvin-Helmholtz instabilities have formed. 



3.3 Self-Gravity 

This test follows, in three dimensions, the collapse of an 
initially uniform sphere to form an n = 1 polytrope. A poly- 
tropic cloud is one in which the pressure P varies according 
to the following relation: 
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Table 5. Ray leigh- Taylor parameters. 



Variable (Unit) 


Value 


Description 


pi (gcm-2) 


2 


Surface density of upper material 


P2 (gcm-2) 


1 


Surface density of lower material 


g (cms-2) 


0.1 


gravitational acceleration 


E.O.S. 


Adiabatic 


Equation of state 


L (cm2) 


1 


Computational domain size 


'^cells 


5122 


Number of grid cells 



Table 6. Self-gravity test parameters. 



Variable (Unit) 



Value 



Description 



-^sphere 
^sphere (pc) 
7 

E.O.S. 

K 

Hcells 



1 
1 

2 

Polytropic 
4.1317 X 1029 
128^ 



Mass of initial sphere 
Radius of initial sphere 
Adiabatic index 
Equation of state 
Equation [l2l constant 
Number of grid cells 



(12) 



where n is the index of the polytrope, p is the density and 
X is a constant. 

The corresponding solution to Poisson's equation for a 
self-gravitating polytropic fluid is given by the Lane-Emden 
equation, which details the variation in pressure and density 
in terms of dimensionless variables C and 0: 



]__d_ 

C2 dC 



0. 



and C are given by equations |14| and |15| 

pc 

and 



C = r 



(n + l)Pc 



(13) 



(14) 



(15) 



where r, pc and Pc are the radial position, central density 
and pressure respectively. 

The solution to the Lane-Emden equation for an n = 1 
polytrope is simply a sine function, which should be the 
form of the resulting density distribution once collapse has 
occurred. We employ reflecting boundary conditions. A sum- 
mary of the parameters used for this test is given in Table |6] 
We ran the model with significant artifical viscosity in order 
to strongly damp the oscillations that would otherwise oc- 
cur. The resulting radial density distribution as calculated 
both analytically and by TORUS is given in Figure |6] and 
demonstrates that TORUS is in excellent agreement with 
the expected result. 



3.4 Radiation Hydrodynamics: Expansion of an 
HII region 

The time evolution of the extent of an HII region around, 
for example, an 0-type star can be divided into two regimes. 
First note that using a time-independent, photoionization 
equilibrium consideration, the HII region extends until the 



CO 1e-22 

'e 



TORUS 
Analytical 



0.2 0.3 0.4 0.5 0.6 

Distance, pc 



Figure 6. The analytical and TORUS-computed resulting den- 
sity distribution for the collapse of a uniform density sphere to 
form an n = 1 polytrope. 



rates of ionization and recombination become equal. In a 
uniform density, hydrogen-only, isotropic medium this cor- 
responds to a spherical HII region with extent known as 
the 'Stromgren radius', given by: 



ri 



( 3iV^ \ 
V47rn2a(2); 



1/3 



(16) 



Where is the number of ionizing photons from the source 
per second, Ue is the electron number density and a^^-* is the 
recombination coefficient into all states except the ground 
state. 

In the time-dependent case the first phase of evolution 
is that in which ri < r^. Here, the expansion is initially 
very rapid. Since the surrounding material near the star can 
be assumed to be ionized as soon as the stellar radiation 
reaches it, the ionization front propagates at the speed of 
light over a mean free path. Given that the ionizing front 
propagates much more rapidly than the speed of sound Cg 
little subsequent material motion occurs. 

In the second phase of evolution, the hot ionized re- 
gion expands into the surrounding, cool, material as a con- 
sequence of the pressure difference between them. Once the 
ionization front expansion velocity drops below Cg, a shock 
moves ahead of the ionization front into the neutral mate- 



rial. Spitzer ( 1978 ) showed that the analytical radius of an 



HII region in the phase two expansion at time t is given by 

4/7 



n 



rn 1 + 



7 cit 
4 rf 



(17) 



up until the point where pressure equilibrium is approached, 
where ci is the speed of sound in the ionized gas. Equation 
|17| is constructed using the thin shell approximation. 

In this test the second phase expansion is modelled in 
three dimensions, using equation ^| as a comparison. The 
evolution of the first phase expansion, prior to r^, is ignored 
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Table 7. Parameters used for the HII expansion model. 



Variable (Unit) 


Value 


Description 


Po (mHcm-3) 


100 


Initial density 


(K) 


10 


Initial temperature of grid 


Uo (cms -*■) 





Initial velocity throughout grid 


7 


5/3 


Adiabatic index 


E.O.S 


Isothermal 


Equation of state 


n (K) 


40000 


Effective source temperature 


(Ro) 


10 


Source radius 


L(pc3) 


11.36 


Grid size 


^cells 


128^ 


Number of grid cells 



tions at a given point in time and has also been shown 
to produce Kelvin- Helmho It z and Ray leigh- Taylor instabili- 
ties. Furthermore our self- gravity calculation reproduces the 
same n = 1 polytropic density distribution as given by the 
Lane-Emden equation following the collapse of a uniform 
density sphere. Finally, the radiation hydrodynamics scheme 
has been shown to agree with the analytical work of |Spitzer| 
([1978) for the rate of expansion of an HII region. These 
results verify that TORUS hydrodynamics and photoioniza- 
tion modules reproduce the standard radiation and hydro- 
dynamical benchmarks and is therefore ready for application 
to new problems. 



2.4 r— 
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1.6 - 
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1.2 - 
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0.8 i 

0.6 r 
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Analytical ■ 



0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 
Tinne, kyr 

Figure 7. The position of the ionization front with time from 

n = r° 



because the calculations here assume photoionization equi- 
librium. The system consists of a star at 40000 K with a 
blackbody spectrum at the centre of a 4pc^ box of neutral 
hydrogen with reflective boundary conditions. We perform 
a radiation hydrodynamics calculation as outlined in section 
|2] and follow the evolution of the ionization front position, 
defined as the point where the atomic hydrogen ionization 
fraction X(HI) = 0.5, from n = with time. 

Table [7| lists the parameters used in this test and the 
results are shown in Figure [7| TORUS shows excellent agree- 
ment with equation [T7| shortly after reaching the Stromgren 
radius. (The discrepancies at early times occur because 
TORUS evolves from a neutral starting point whereas the 
analytical solution starts with the ionization front at the 
Stromgren radius). At late times the evolution starts to de- 
viate from the analytical solution as sufficient material is 
accumulated for the thin shell approximation to no longer 
apply. 



3.5 Testing Summary 

We have shown that TORUS satisfies a number of bench- 
mark tests. The radiative transfer scheme, including treat- 
ment of the diffuse field, is in good agreement with the Lex- 



ington benchmark as calculated by Cloudy (Ferland et al 



4 RADIATIVELY DRIVEN IMPLOSION OF A 
BONNOR-EBERT SPHERE 

Here we model the radiatively driven implosion of a Bonnor- 
Ebert sphere (a sphere in which the density varies according 



to equation 13) using three different treatments of the ra- 



diation field so as to distinguish their relative contributions 
to the evolution of the system. The three different radiation 
fields used are: 

a) a monochromatic radiation field with the OTS ap- 
proximation 

b) a polychromatic radiation field with the OTS ap- 
proximation 

c) a polychromatic radiation field with the diffuse field, 
as outlined in section |231 

The details of the system are very similar to that of 
Gritschneder et aTr|p()09| ). A Bonnor-Ebert sphere of radius 



1.6 pc resides at the centre of a 1.5 x 10 cm (4.87 pc ) grid 
of 128^ cells. Our model domain is slightly larger than that 
used by Gritschneder et al. (12009) so that the evolution of 
larger extent of any shadowed region can be studied. The 
BES has a core number density of lO^cm"^, with the ma- 
terial surrounding the BES having a number density equal 
to that at the BES edge (such a BES is known as 'incom- 
plete'). The resolution of these models is currently limited 
by computational cost, however the number of grid cells here 
is equivalent to that used in the successful Hll-region expan- 
sion test of section [3^ This model is also on a smaller length 
scale than that of the Hll-region expansion test so the reso- 
lution will be sufficiently high. As mentioned in section [2^ 
we will be implementing an AMR grid that will enable the 
calculation of models at higher resolutions in future work. 

The primary photon source is a star that lies outside the 
grid in the —x direction. The radiation field is assumed to 
enter the grid plane parallel at the —x boundary, at which 
photon packets are initiated at random locations and the 
flux is modified to account for geometric dilution. 

As well as considering the three different radiation 
schemes mentioned above, we also treat the three different 



flux regimes considered in Gritschneder et al. (2009). These 



1998). The hydrodynamics algorithm satisfies the Sod shock 



tube test and Sedov- Taylor blast wave density distribu- 



are denoted high, medium and low flux and correspond to 
the BES being located just within, on the edge of and just 
beyond the Stromgren radius respectively. The fluxes and 
corresponding stellar properties that were used to generate 
these different flux regimes are given in Table [S] along with 
the other parameters used for this model. 

In all of the models presented here, the hydrody- 
namic boundary conditions are periodic at ±y and ±z 
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Table 8. Parameters used for the RDI of a BES model. 



Variable (Unitj 


Value 


Description 


Rc (pc) 


1.6 


Cutoff radius 


^max (cm~^) 


1000 


Peak BES number density 


To (K) 


10 


Initial temperature of grid 


7 


1 


Adiabatic index 


E.O.S 


Isothermal 


Equation of state 


(cm-2) 


9.0 X 10^ 


Low ionizing flux 


Ao (pc) 


(-10.679,0,0) 


Source position (low flux) 


^rned (cm-^) 


4.5 X 10^ 


Intermediate ionizing flux 


^med (pc) 


(-4.782,0,0) 


Source position (medium flux) 


$hi (cm-2) 


9.0 X 10^ 


High ionizing flux 


^hi (pc) 


(-3.377,0,0) 


Source position (high flux) 


L(pc3) 


4.87 


Grid size 


''^cells 


128^ 


Number of grid cells 


CFL 


0.3 


CFL parameter 



do not include periodic photon packet boundary conditions 
and so these regions are subject to a non-symmetric diffuse 
ionizing flux. 

The logarithmic density distribution for the high, 
medium and low flux models over all three treatments of 
the radiation field, are shown side by side at 50, 100, 150 
and 200 kyr in Figures |9] [lO] and [TT] respectively. In each 
of these figures the monochromatic models are represented 
by the left hand column, the polychromatic models by the 
central column and the polychromatic-diffuse models by the 
right hand column. It is clear that there are some marked 
differences between the evolutions of the system under the 
different radiation treatments. A case-by-case study of the 
evolution of the separate models is given below in sections 
|4T3l[4T4l and [4T5l 



and free outflow/no inflow at zbx. Dirichlet boundary 
conditions are used for the self-gravity calculation and 
the radiation field boundaries are free outflow/no inflow. 
Using this boundary condition for the radiation field 
can lead to reduced sampling at the domain boundaries, 
where material will be subject to non-symmetric diffuse flux. 

The free fall time for this cloud is approximately 3 Myr, 
estimated using l/(^G'pmax) where pmax is the central den- 
sity and G is the gravitational constant. This is a factor of 
15 longer than the total simulation time of 200 kyr. Radia- 
tion hydrodynamics will thus dominate gravitational effects 
in the evolution of the system. We dump the state of the 
computational grid every 5 kyr. 



4.1 Results and Discussion 

4-1- 1 Initial properties 

The initial ionization state of the system for all three flux 
regimes and photoionization schemes is shown in Figure 
[S] At this stage there is already a noticeable difference 
between them in the extent of their un-ionized regions. The 
top row from Figure [S] represents the starting point that 
would be obtained by most pre-existing models (specifically, 
it is a zoomed out version of Gritschneder et al. 2009). 



The middle row is the start point for the models in which 
a polychromatic radiation field is considered, but the 
GTS approximation is still applied. In comparison with 
the top row, it is clear that the extent of the ionized 
region has increased slightly and the transition region 
has been smoothed out. This is due to hard radiation, 
which penetrates more deeply since the photoionization 
cross section approximately decreases in proportion to the 



inverse cube of the photon frequency (Osterbrock 1989) 



The bottom row is the starting Hi fraction for the models 
that include the diffuse field. In all three flux regimes the 
extent of the ionized region is significantly different to that 
of the other two sets of models. At high flux material in the 
wings of the model is completely ionized, the medium flux 
model I-front has significantly wrapped itself around the 
BES and the low flux model I-front now grazes the BES. 
Note that the curved I-front wings towards the edge of the 
low and medium flux models arise because these models 



4.1.2 The formation and evolution of instabilities 

A number of the models presented here are subject to the 
formation of instabilities. These arise when the thin shell 
swept up by the ionization front propagates though the lower 
density regions of the computational domain such as in the 
wings or (in the low flux regime) prior to driving into the 
BES. 

There are two main sources of perturbation that may 
seed these instabilities. The first is 'angle of incidence' 
Williams ( 2002 ) in which regions of the I-front that are not 



entirely perpendicular to the incoming radiation field are 
subject to varying ionizing flux and therefore differential 
acceleration. This may occur as the I-front wraps around 
the BES. The second is numerically induced perturbation 
via noise in the calculated ionization fraction. As mentioned 
in section |2.3| the induced temperature, and therefore 
pressure, is directly proportional to the calculated ioniza- 
tion fraction (equation [7| . The resulting pressure gradient 
then determines the induced advecting velocity. If a thin 
shell of material does not encounter any disruption in its 
propagation (like encountering a high density component 
of a BES) then even a small amount of numerical spread 
in the ionization fraction along the I-front will eventually 
lead to it bending on small scales and therefore induce 
thin shell instabilities (Vishniac 1983| Garcia-Segura & 
Franco 1996). Once the I-front structure is disrupted, the 



faster propagating components will lose mass by transfer 
to the cool neutral neighboring material perpendicular to 
the direction of the I-front propagation. This leads to an 
accumulation of material ahead of the slower moving com- 
ponents of the I-front, which further brakes the expansion 
in these regions. This transport of material also reduces 
the density in the faster moving components of the I-front, 
allowing photoionizing radiation to propagate more deeply 



(c.f. equation 16) and accelerate these components of the 
front further. Improving the accuracy of the ionization 
fractions in the front will delay the onset of numerically 
induced instabilities, however they should eventually arise 
for any non-analytical radiation hydrodynamics code if the 
I-front is allowed to propagate for long enough without 
being disrupted by some other means. 

It is clear that other radiation-hydrodynamical methods 
should also seed these instabilities. For example, combining 
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Figure 8. Top row) The hydrogen ionization fraction for the OTS monochromatic model. Middle row) The hydrogen ionization fraction 
for the OTS polychromatic model. Bottom row) The hydrogen ionization fraction for the model which includes the diffuse field. Columns 
are, from left to right, high, medium and low flux regimes. Each frame is a slice through the computational grid, which is a cube with 
sides 4.87 pc long. Major ticks are separated by Ipc. 



SPH hydrodynamics with ray-tracing will lead to a 'noisy' 
I-front due to the random-variation of the SPH representa- 
tion of the density field. Grid-based codes with ray-tracing 
radiation-transfer will also be susceptible to instabilities as 
the angular sampling of the radiation field may not coin- 
cide perfectly with the axes of the underlying hydrodynam- 
ical grid. Of course within star forming regions themselves 
density perturbations will inevitably lead to the growth of 
instabilities. 

Regardless of the seed of these instabilities in the sim- 
ulations, their evolution occurs in a manner consistent with 
the instability studies referenced above and result in 'ele- 
phant trunk' structures. Their evolution also highlights in- 



teresting differences between the different treatments of the 
radiation field, the details of which are discussed in the fol- 
lowing sections. 



4- 1.3 High flux models 

The high flux density evolutions are given in Figure |9] and 
broadly exhibit two different behaviours. In the polychro- 
matic OTS and polychromatic-diffuse models the system is 
initially ionized to the extent that a strong shock cannot 
form quickly enough to effectively drive into the BES. What 
material is accumulated propagates for a short time until the 
higher density BES brakes it and a photo-evaporative flow is 




Figure 9. The high flux model logarithmic density distributions (cgs). Left column) Monochromatic models. Middle column) Polychro- 
matic models. Right column) Polychromatic-diffuse models. Time is increasing from top to bottom, with snapshots at 50, 100, 150 and 
200 kyr. Each frame is a slice through the computational grid, which is a cube with sides 4.87 pc long. Major ticks are separated by 1 pc. 
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Figure 10. The medium flux model logarithmic density distributions (cgs). Left column) Monochromatic models. Middle column) 
Polychromatic models. Right column) Polychromatic-diffuse models. Time is increasing from top to bottom, with snapshots at 50, 100, 
150 and 200 kyr. Each frame is a slice through the computational grid, which is a cube with sides 4.87 pc long. Major ticks are separated 
by 1 pc. 




Figure 11. The low flux model logarithmic density distributions (cgs). Left column) Monochromatic models. Middle column) Polychro- 
matic models. Right column) Polychromatic-diffuse models. Time is increasing from top to bottom, with snapshots at 50, 100, 150 and 
200 kyr. Each frame is a slice through the computational grid, which is a cube with sides 4.87 pc long. Major ticks are separated by 1 pc. 
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established ( [Bertoldi McKee||1990| ). The evolution of the 
resulting structure is then a consequence of rocket motion as 
heated, dense, material is evaporated away from the surface 
of the cloud into the low density external material (Oort 



& Spitzer||1955). The tunnelling of material at the tip and 



along the length of these models into the cometary structure 
occurs where there are differences between rocket-motion ve- 
locities due to either variations in the accumulated density 
or the density internal to the shell. An interesting difference 
between the OTS and diffuse field models is also revealed in 
the rocket-driven phase, with the OTS model being acceler- 
ated only along components facing the ionizing source and 
the diffuse field model being accelerated across the entire 
cometary surface due to diffuse- driven photo-evaporation. 

On the other hand, the monochromatic OTS model does 
form a strong shock sufficiently rapidly to effectively drive 
into the BES, this leads to greater compression and accumu- 
lation of material into a relatively thick shell around the edge 
of the resulting bow structure. To illustrate the early braking 
of the polychromatic models, the difference in the velocity 
field between the monochromatic OTS and diffuse models at 
50 kyr is illustrated in Figure[T2] At this point the monochro- 
matic model is still driving a shock into the BES and accu- 
mulating material, whereas the polychromatic-diffuse model 
is beginning a photo-evaporative flow. 

In the OTS model, suflicient material is rapidly accu- 
mulated for a short, strong, photoevaporative flow to oc- 
cur prior to substantial braking of the bow. The resulting 
rocket-motion is therefore much stronger than normal photo- 
evaporative flow, giving rise to the ejection of a signiflcant 
amount of material that accelerates the existing shock and 
carves out a low density wake. Subsequent ejections also oc- 
cur episodically along the length of the bow that is exposed 
to ionizing radiation. The result is a disrupted region sur- 
rounding the tip of the bow structure in which densities can 
be excavated to levels lower than the ambient surroundings. 
A possible cause of the episodic nature of this process is 
that the outer shell density oscillates about some critical 
value as the shell sequentially accumulates and ejects. This 
'episodic photo-evaporative ejection' further drives the col- 
lapse very effectively and contracts the bow perpendicularly 
to the ejection direction, tapering the head of the cometary 
structure. Figure [13] shows the disrupted region around the 
tip of the bow of the monochromatic model at 180, 185 and 
190 kyr and illustrates the motion of discrete knots of mate- 
rial away from the surface, rather than a continuous stream. 

The evolution of the maximum cell densities for each 
high flux model are shown in the top frame of Figure |14| 
This can be used in conjunction with the appropriate density 
map to illustrate the rate at which material is collected. Here 
the maximum density evolution highlights the differences be- 
tween the two different behaviours noted above. In the poly- 
chromatic and polychromatic-diffuse models the maximum 
density increases at a declining rate and eventually plateaus. 
The monochromatic model continues accumulating material 
as it is effectively rocket-driven towards the core of the BES, 
flnishing with a maximum density approximately 4.5 times 
that of the other models. Star formation at this flux regime 
will actually occur more slowly when polychromatic radia- 




Figure 12. The logarithmic density field in gcm""^, with velocity 
vectors, for the high flux models at 50 kyr. This demonstrates 
the difference between the driving shock of the monochromatic 
OTS model (top frame) and the photo-evaporative flow of the 
polychromatic-diffuse model (bottom frame). Each frame is a slice 
through the computational grid, which is a cube with sides 4.87 pc 
long. Major ticks are separated by 1 pc. Typical shock (upper 
frame) and outflow (lower frame) velocities are 5-7 kms"-*^ and 
1-4 km s"-*^ respectively. 



tion or the diffuse fleld is accounted for than in the simplifled 
calculation. 

The formation of thin shell instabilities has little to no 
impact in the evolution of these models, appearing only in 
the wings of the polychromatic and monochromatic OTS 
models and are simply swept off the grid. 



4- 1-4 Medium flux models 

The medium flux density evolution is presented in Figure 
and exhibits the largest difference between the flnal states 
of the OTS and diffuse fleld models. The OTS models both 
develop shocks that drive effectively into the BES, resulting 
in rapidly formed high density bow structures. The accumu- 
lated density is high enough to give rise to a scaled down 
version of the episodic photo-evaporatively driven collapse 
exhibited by the high flux monochromatic model and leads 
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Figure 13. Logarithmic density snapshots in gcm""^ with veloc- 
ity vectors for the high flux monochromatic model at 180, 185 and 
190 kyr. This illustrates the ejection of distinct knots of material 
and the evolution of the disrupted region around the tip of the 
bow. Each frame is a slice through the computational grid, which 
is a cube with sides 4.87 pc long. Major ticks are separated by 
1 pc. Typical velocities in the outflow region of each frame are 
25-35 km s-^ 
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Figure 14. The evolution of maximum logarithmic densities re- 
spectively. Flux regimes are high to low from top to bottom. 



to the same rocket-motion that continues driving material 
towards the centre of the BES and tapers the bow head. 
In the wings of the OTS models instabilities form via the 
mechanism described in section [4.1.2| and propagate linearly, 
having no effect of the rate of collapse and eventually being 
swept off of the grid. Note that the low density regions in the 
wake of these instabilities are not due to photo-evaporative 
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flow, rather they have been excavated by the high velocity 
instability shocks, with the material funnelled laterally to 
form knots and elephant trunks. 

The diffuse field model behaves slightly differently, 
resembling a more effectively collapsed version of the 
monochromatic OTS high flux model. Again a strong shock 
is developed which drives into the BES and transitions to 
episodic photo-evaporatively driven collapse. The major dif- 
ference lies in the wings of the model, which effectively drive 
into the side of the BES resulting in a final bullet-like struc- 
ture with a very high density at what was the core of the 
BES. The reason for this is that a photo-evaporative flow 
can establish itself along a greater extent of the bow, into 
regions that would otherwise be shadowed from heating, be- 
cause of the diffuse radiation field. This results in a whip-like 
progression of photo-evaporation along the trunk as the shell 
sequentially becomes dense enough to readily eject material 
and ceases in the tail regions when the diffuse ionizing flux 
becomes too low to cause heating. This will most likely re- 
sult in much more rapid star formation in the model that 
includes the diffuse field. The evolution of the maximum cell 
density for these medium flux models is shown in the middle 
frame of Figure This clearly highlights the fairly extreme 
difference between the OTS and diffuse field models, with a 
final difference in maximum density at 200 kyr of well over 
a factor of 10. 



4-1-5 Low flux models 

The low flux density evolutions are given in Figure[lT] These 
models exhibit the largest susceptibility to thin shell insta- 
bility as instabilities can begin developing across the en- 
tirety of the I- front before it impacts the BES. Instabilities 
that collide with the BES are braked and result in high 
density knots along the rim of the resulting bow structure. 
Those that continue in the wings of the model are elon- 
gated into elephant trunks with high density tips via the 
mechanism described in sect ion 14.1.2 1 The evolution of these 
elephant trunks varies with each treatment of the radia- 
tion field. The polychromatic OTS trunks are more elon- 
gated than the monochromatic ones because hard radia- 
tion carves out a path more rapidly. Any compression of 
the trunks in the OTS models is due to thermal pressure. In 
the polychromatic-diffuse model, diffuse field radiation effec- 
tively drives into the material perpendicular to any displace- 
ment in the ionization front and therefore actually prevents 
the formation of a number of potential trunks by smoothing 
out dimples. Those trunks that do form will also continue to 
be both compressed thermally and exposed to diffuse ioniz- 
ing radiation. 

The RDI process for the OTS models occurs very 
weakly, with more material being accumulated through in- 
stability than compression. The diffuse field model drives 
into the BES more effectively, forming a smoother high den- 
sity bow compared to the knotted structures of the other 
models and manages to initiate some photo-evaporative 
flow. The maximum cell density evolution for the low flux 
models is given in the bottom frame of Figure At 200 kyr 
the diffuse field model has accumulated only a slightly larger 
maximum density than the OTS models, though it is clear 
from Figure^^that it has achieved this order of density over 
a relatively large volume. The diffuse field model will most 



likely form stars first and on a larger scale than the models 
that use the OTS approximation. 

These low fiux results, comprising a bright rimmed 
cloud with pillars along its wings, bear resemblance to ob- 
served systems such as IC 1848E, as shown in Figure 1 of 



Chauhan et al. (2011). In the aforementioned work, insta- 



bility was suggested as the formation mechanism of the ele- 
phant trunks in this region of IC 1848E, our unstable radia- 
tively driven shock driving into a pre-existing density struc- 
ture supports this hypothesis. 



4-1-6 Comparison with iVINE 

There are some differences between the results obtained here 
and those of [Gritschneder e t al.| (j2009|, the work on which 
our model parameters are based. In particular the maximum 
density evolutions derived by TORUS are much weaker, to 
the extent that no clear gravitational collapse has occurred 
by 200 kyr. This discrepancy is attributed to the difference 
in grid size, combined with the use of periodic boundary 
conditions. A significant proportion of the compression of 
the BES in the iVINE models arises because hot gas is ad- 
vected off the edge of the periodic boundary and impacts 
the cloud on the opposite side of the grid (see section 4.1.2 
of Gritschneder et al. 2009). This effect is justified by the 
authors as they assume that the molecular cloud completely 
surrounds the triggering star, however the effect will still 
give rise to differences in the results obtained by our models 
and theirs since motion of the hot gas has longer to decay 
over our larger grid and is also impeded by the elephant 
trunks that have formed in the wings of our models. In our 
models, lateral compression of the BES due to motion in the 
hot gas through the boundaries has negligible effect and as 
such the BES can be considered to be isolated. 

With regard to the effect of the diffuse field, a compar- 
ison with the results found using iVINE/DiVINE ( [Ercolano" 
&: Gritschneder|2011b ) is not straightforward as the systems 
being modelled are very different. However we broadly agree 
that treating the diffuse field can lead to higher density re- 
sulting structures, particularly in our medium fiux model 



(section 4.1.4). In our low fiux model (section 4.1.5), where 
elephant trunks of a similar form to those created in the 
iVINE/DiVINE models arise, we also find that inclusion of 
the diffuse field increases the density and decreases the num- 
ber of elephant trunks that form. We also agree that these 
trunks are narrowed, sometimes to the extent that the head 
of the trunk can become completely detached. 



5 CONCLUSIONS 

It is clear that a direct treatment of the diffuse field has 
a significant impact upon the evolution of radiation hydro- 
dynamics calculations on parsec scales. Not only is the ef- 
fectiveness of RDI sensitive to the way in which the radia- 
tion field is treated, the formation and evolution of elephant 
trunk structures via instability also varies. 

At low and intermediate fiux regimes inclusion of the 
diffuse field results in more efficient RDI, with both more 
widespread accumulation of material and, particularly in 
the medium fiux case, a higher final maximum density. In 
the high fiux regime however, the material around the BES 
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is ionized so rapidly in the diffuse and polychromatic OTS 
models that only a weak shock forms and RDI does not 
occur effectively. Generally, it is found that the extent to 
which RDI occurs depends strongly on the strength of the 
shock that is accumulated prior to driving into the BES. 
Regardless of the ionizing flux, if only a weak driving shock 
has formed collapse will occur slowly. When sufficient ma- 
terial is accumulated a photo-evaporative flow is found to 
occur before the driving shock is braked, reinforcing the 
shock with the resulting rocket motion. This occurs in short, 
sharp, bursts with the reason for this episodic behaviour hy- 
pothesised as being due to the refilling time for material in 
the dense shell. The high velocity ejecta carve out material 
around the head of the cometary structure resulting in a low 
density disrupted bubble. This rocket motion has a signifi- 
cant effect on the collapse of the BES: with the inclusion of 
the diffuse field, it can accur across a large part of the other- 
wise shadowed region and significantly narrow the resulting 
cometary structure. 

Elephant trunks that arise due to thin shell instability 
are harder to form in the presence of the diffuse field as the 
dimples that seed them are smoothed out, and those that do 
form are then subject to the expected combination of ther- 
mal compression and diffuse field photoionizing radiation. 
Despite this the formation of elephant trunks, as the result 
of instability in a radiatively driven shock, still occurs and 
provides a mechanism for the formation of systems such as 
that discussed in Chauhan et al. (2011). 



The radiation-hydrodynamical effect of the diffuse 
field is cumulative and significant, and is strongly coupled 
to the hydrodynamical evolution of the system. This 
implies that the treatment of the diffuse field should be 
accurate throughout a simulation, as consistent deficiencies 
could lead to a systematic change in the overall radiation 
hydrodynamical evolution of the system. 

We next intend to compute synthetic images to deter- 
mine theoretical observables for systems undergoing star for- 
mation as a result of radiative feedback. We will also include 
the use of an adaptive mesh so that we can follow collapse 
until the onset of star formation. 



6 SUMMARY 

We have used the radiation hydrodynamics module of the 
TORUS code to investigate the effects of using a monochro- 
matic OTS, polychromatic OTS and polychromatic-diffuse 
field on the radiatively driven implosion of a Bonnor-Ebert 
sphere. We have found that incorporating the diffuse field 
into this model over three flux regimes leads to significantly 
different results to those obtained using the OTS approx- 
imation. At intermediate and low flux regimes the rate of 
compression is higher than that without inclusion of the dif- 
fuse field, whereas at high flux compression actually occurs 
more slowly when the diffuse or polychromatic OTS field is 
considered because there is little opportunity for a material 
shock to drive into the BES. In the event of accumulation 
of sufficient material, photo-evaporative flow or ejection has 
been identified as a mechanism which can drive collapse very 
effectively. This photo-evaporative flow is particularly effec- 
tive at driving and compressing the tail of the cometary 



structure when the diffuse field is treated. The formation of 
elephant trunk structures via instability also occurs much 
less readily with the inclusion of the diffuse field as pertur- 
bations to the ionization front are smeared out. We conclude 
that in order to properly address quantitative questions re- 
garding triggered star formation thorough treatment of the 
diffuse field is necessary in radiation hydrodynamics models. 
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